Thermodynamics of C incorporation on Si(100) from ab initio calculations 
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We study the thermodynamics of C incorporation on Si(100), a system where strain and chemical 
effects are both important. Our analysis is based on first- principles atomistic calculations to obtain 
the important lowest energy structures, and a classical effective Hamiltonian which is employed to 
represent the long-range strain effects and incorporate the thermodynamic aspects. We determine 
the equilibrium phase diagram in temperature and C chemical potential, which allows us to pre- 
dict the mesoscopic structure of the system that should be observed under experimentally relevant 
conditions. 
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PACS numbers: 61.66.Dk, 68.35.Rh, 68.35.Bs 

Carbon-enriched silicon systems are the focus of cur- 
rent interest as candidates for a material with tailored 
electronic properties, which is compatible with well- 
established silicon technology. The tetravalent nature 
and large band gap in the diamond structure make car- 
bon an ideal candidate for incorporation in Si. However, 
the solubility of C in Si under thermodynamic equilib- 
rium is extremely low (« 10 -5 ) due to the huge mismatch 
in bond length (35%) and bond energy (60%) between 
C and Si. Non-equilibrium methods, such as molecular 
beam epitaxy (MBE), that exploit the higher atomic mo- 
bility on surfaces, can be used to overcome this obstacle 
and enhance solubility [Q . As predicted theoretically by 
Tersoff^], C solubility is enhanced by several orders of 
magnitude near the Si(100) surface, especially in sub- 
surface layers. Osten et al. || confirmed experimentally 
this prediction and observed that C atoms diffuse to sub- 
surface layers above a certain temperature. This finding 
opened new possibilities for growth of C-rich metastable 
structures. 

The enhanced solubility near the surface has impor- 
tant consequences. The large tensile strain associated 
with C incorporation in Si has proven a very powerful 
tool in device engineering: a small amount of C can com- 
pensate the Ge-induced tension in pseudomorphic SiGe 
layers H , and can also suppress dopant outdif- 

fusion. This idea was recently implemented in a novel 
heterojunction bipolar transistor (?| . Another interesting 
effect^ [| produced by C incorporation on the Si(fOO) 
surface is an unusual change of the surface periodicity 
after deposition of even a small amount of C («s |of a 
mono-layer (ML)): the well-known c(2 x 4) orp(2 x 2) re- 
constructions of the pure Si surface change to a c(4 x 4) 
pattern. This is clearly visible in several LEED experi- 
ments after ethylene exposure Jf(], [ll[] , or MBE j|] . 

The microscopic features of C incorporation in Si are 
rather well understood. Previous work by the authors 
Jf2| , [fg| , |l4| revealed an oscillatory C profile driven by 



the competition between two factors: the tendency of 
C atoms to occupy favorable sites which are determined 
by the reconstruction strain field, and the preferential 
arrangement of C atoms at certain distances which min- 
imizes the lattice elastic energy. The profile is character- 
ized by enhancement of C content in the first and third 
layer, depletion in the second and an exponential reduc- 
tion from the fourth layer and beyond. The favorable C 
sites are in the third layer between the surface dimers, be- 
ing under compressive stress ^ and thus suitable for 
the smaller-sized C atom. The preferential arrangement 
of two C atoms in the surface layers is at a third nearest- 
neighbor position, which is also the lowest-energy config- 
uration in the bulk Jf| . The interaction of C atoms in first 
nearest neighbor position is highly repulsive |]l2| . More- 
over, a recent experiment || showed that C-C dimers are 
quite rare in the annealed C-Si(lOO) surface. 

On the other hand, a mesoscopic picture of the surface 
structure that would link the atomistic features to long- 
range strain effects, and predict structures for different 
growth conditions more relevant to experiment, is lack- 
ing. Here, we present an approach that closes this gap 
and is able to determine the equilibrium surface phase 
diagram. It is based on first-principles atomistic stud- 
ies, incorporating microscopic strain and polarity effects, 
which are expected to have a dominant role in the under- 
standing of C incorporation in the host Si lattice. The 
results are linked to a classical effective Hamiltonian in 
a Monte Carlo scheme that incorporates strain into the 
thermodynamic aspects of the problem. 

We considered first all possible atomistic configura- 
tions likely to have low energies. Guided by previous 
work, ^ [f2], [f3| [f(| we can establish a set of rules which 
should be obeyed by low energy configurations of the C- 
Si(IOO) system, namely: 

(i) The surface has c(4 x 4) periodicity. 

(ii) C substitutes Si atoms only in the first (surface) or 
in the third layer. 
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FIG. 1: Top views of the relaxed geometries of the 12 c(4 x 4) 
cells we considered. Si surface atoms are shown as grey circles, 
Si atoms in 2nd and 3rd layers as smaller, white circles, and 
C atoms as black circles. 



(iii) There are no C-C nearest neighbors, or second 
nearest neighbors. 

(iv) Third layer C exists only at sites between pairs of 
surface dimers. 

We constructed every possible structure of the C-Si(100) 
system, consistent with these rules, with 0, 1, 2, 3 and 4 
C atoms per c(4 x 4) unit cell. These are shown in Fig. 1. 
The configurations are named nX, n being the number 
of C atoms in the c(4 x 4) unit cell and X an index to 
distinguish structures of the same n. 

To study the energetics and the geometrical fea- 
tures of these configurations, we performed ab initio 
DFT/LDA/Pseudopotential calculations 0. In Table I, 
we give the energy E(nX) for each configuration nX, rel- 
ative to the pure Si(100) surface (configuration 0). These 
energies are defined by: 

E{nX) = E tot {nX) - n(^ c - Hsi) - E tot (0), (1) 

where E to t(nX) is the calculated total energy of config- 
uration nX, and /icsMSi are the chemical potentials for 
C and Si atoms. For Si, the chemical potential is taken 
to be the energy of the bulk, since this is the natural 
reservoir for Si atoms due to the presence of steps, ter- 
races, and other surface defects. As a starting point, [ic 
is taken to be the energy of a C atom in bulk diamond. 
We discuss its variation below. 



TABLE I: Relative energies for the 12 configurations of Fig. 
1, according to Eq. (1) in eV per c(4 x 4) unit cell, nc is 
taken as the energy of a C atom in diamond. ]l8[ 



For low concentration of C, one atom per unit cell, 
which corresponds to | ML coverage, the structure with 
subsurface C (IB) has lower energy than the one with a 
Si-C dimer (1A) on the surface. This result is in agree- 
ment with experiment: X-ray spectra |l(| indicate that 
for low C deposition, there is almost no C on the sur- 
face. A single C atom at a third-layer site relieves the 
compression due to the surface reconstruction and makes 
four Si-C bonds while a single C atom on the surface, as 
part of a dimer, allows the formation of only three SiC 
bonds. Passing to configurations with more than one 
subsurface C atom per unit cell, we find that the energy 
of configuration 2B is considerably higher than that of 
2A, 2C, and 2E. Apparently, the presence of two sub- 
surface C atoms per unit cell produces a large distortion 
of the surrounding Si-Si bonds. Thus, for two C atoms 
in the unit cell (4 ML coverage), the best situation is to 
have both C atoms on the surface. For three C atoms per 
unit cell (| ML coverage), we observe an equivalence of 
surface/subsurface sites: configurations 3 A and 3B have 
almost the same energy. For the case of four C atoms ( \ 
ML coverage), there is just one structure (4A) consistent 
with the rules we discussed before. For illustration pur- 
poses, we include another configuration, 4B, containing 
C atoms in second neighbor positions, which violates rule 
(iii); its very high energy can be considered as a justifica- 
tion for these rules. We do not go beyond \ ML coverage 
because experiments show that this results in disordered 
structures. 

To compare configurations with different C content, we 
use Eq. (1). This implies that the surface is in equilib- 
rium with a reservoir of C atoms characterized by \ic- It 
is generally accepted that C forms small clusters, whose 
cohesive energies are in the range between —5.5 and —7 
cV [|l9|. In Fig. 2(a), we show the energies of configu- 
rations 0, IB, 2A, 3A and AA as a function of fic- The 
vertical faint lines denote the transition points at which 
the preferred structure changes. We observe that the low- 
est energy structure depends strongly on [ic and hence 
on the conditions of C deposition, making it difficult to 
predict what the actual structure of the system will be. 
This implies that the equilibrium surface structure might 
be composed of different configurations at the atomistic 
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FIG. 2: (a) The energy of the preferred structure for 1,2,3 
and 4 C atoms per unit cell as a function of uc, from Eq. 
(1). (b) The average C content of the system as a function 
of C chemical potential, at 850K, with (filled symbols) and 
without (open symbols) the elastic strain interactions. The 
inset shows the ratio of surface to subsurface C atoms at the 
same temperature, as a function of the total C content in the 
unit cell, (points marked with bigger symbols correspond to 
the cases shown in Fig. 3). 



level and, consequently, a larger-scale mesoscopic pic- 
ture is needed. In this picture, the elastic interactions 
in the boundaries between regions of different C content 
or atomistic structure will play an important role. 

To obtain such a mesoscopic description and study its 
thermodynamics, we use an effective Hamiltonian within 
a Monte Carlo scheme. We first construct a classical, 
12-state generalized Potts model. We assume that the 
surface consists of c(4 x 4) cells, each of which can exist in 
any of the lowest energy configurations discussed before. 
The Hamiltonian is 



H 



eff 



(2) 



where the summation on i runs over all N cells, and 
the summation on j over the 8 nearest neighbors of i 
(the factor of one half takes into account double-counting 
of the interaction); Ci is the configuration in cell i, and 
AE(ci,Cj) is the interaction energy between the neigh- 



boring cells Ci and Cj (AE(ci,a) = 0). In order to ob- 
tain these interaction energies, we calculated the total 
energies of all 144 c(12 x 12) cells consisting of a cell 
Ci surrounded by 8 Cj cells. These were obtained using 
Tersoff's empirical potential ]20| ], with a Monte Carlo re- 
laxation method. The key point here is that we are inter- 
ested in elastic strain interactions between different cells 
rather than any specific features of their structure; for 
this type of interaction between cells the classical poten- 
tial employed gives very reasonable results as established 
by previous studies of similar systems || ^3|. Having 
determined 7Y e //, we performed an equilibrium Monte 
Carlo study for different values of the temperature and 
fic- We used a 70 x 70 grid of cells, which was suffi- 
cient to get converged averages. For each value of fic the 
system was started from a random configuration at high 
temperature (1200K) and then was cooled down to room 
temperature. 

In Fig. 2(b) we plot the average total C content of the 
system as a function of fic- The zero for fic is taken to 
be the energy of an isolated C atom. For comparison, wc 
plot the same quantity when the elastic strain interac- 
tions, AE(a,Cj), are set equal to zero. In both cases the 
average C content of the system increases with increasing 
fic- The difference between the two curves shows the im- 
portance of elastic interactions, which, being of the order 
of 0.1-0.2 eV/ c(4 x 4) cell, are more important than the 
total energy differences in the low fic region. This effect 
weakens as fic increases. 

Another important consideration, with experimentally 
observable consequences, is the relative amount of sur- 
face versus subsurface C atoms as a function of the to- 
tal C content Ctot in the unit cell, shown in the inset 
of Fig. 2(b). For larger total C coverage the ratio in- 
creases monotonically, approaching infinity as Ctot — * 4. 
Again, we give the case with AE(ci,Cj) — for com- 
parison. There is a critical coverage of 1.39 C atoms per 
unit cell, or 0.17 ML, beyond which C prefers mostly sur- 
face sites. This critical coverage has a weak temperature 
dependence: it has a maximum around 850K, and then 
falls monotonically down to roughly 0.16 ML at 300K and 
HOOK. The increase of the surface C with increasing to- 
tal C is not unexpected: although third layer compressed 
sites seem ideal for the smaller C atom, a large amount 
of subsurface C would cause large strain and raise the 
total energy dramatically, as exemplified by the case of 
configuration 4i3, which we discussed earlier. From Fig. 
2(b) we see that the transition from mostly subsurface 
to mostly surface C content does not coincide with the 
change from 1 to 2 C atoms per cell, as suggested by 
our ab initio results, indicating that the elastic strain in- 
teractions play a dominant role in the structure of the 
surface. 

Representative snapshots of the C distribution from 
the Monte Carlo run at 850K are shown in Fig. 3. The 
three cases correspond to the points shown in Fig. 2(a). 
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An impressive long range order is revealed for the sub- 
surface C (right column of Fig. 3). For low fj,c, Fig. 
3(c), this order appears in the form of alternating rows 
of similar C content. These consist mostly of unit cells 
corresponding to and 2B structures. Although struc- 
ture 2B has relatively high energy, as shown in Table 
I, a structure of alternating and 2B cells seems to be 
preferable, giving C the chance to take advantage of the 
desirable third layer compressed sites. The line-pattern 
is a result of the strong interactions in that range of /ic, 
which forces the c(4 x 4) cells to have the maximum num- 
ber of different neighbors, 6 for the row pattern and 4 for 
a chess pattern. The latter pattern is observed for higher 
/ic, Fig. 3(b), and is typical for the range of fie and tem- 
perature corresponding to 1-2 C atoms per unit cell on 
average (see also Fig 2(b)). In this case, the self-energies 
of the different configurations dominate over the elastic 
strain interactions. In the specific example shown in Fig. 
3(b), the structure consists mostly of 3A and 2B cells, 
at a ratio 2:1. This suggests that for conditions corre- 
sponding to this range of He, an alternating pattern of 
surface/ subsurface C is preferable. As fic increases the 
energies of the cells become large enough so that the elas- 
tic interactions are not that important, and thus the sys- 
tem does not have much to gain by self-organizing. This 
results in the random pattern shown in Fig. 3(a), which 
consists mainly of 4A cells with a few 3B cells around 
them. Again, of the two configurations with three C 
atoms, the one with some subsurface C is preferred when 
the dominant configuration has all its C atoms on the 
surface. 



Finally, we note that while subsurface C shows this 
interesting ordering behavior, the surface C seems to be 
randomly distributed as the patterns of the middle col- 
umn of Fig. 3 indicate. The ordered structure appearing 
in the total distribution of C (left column of Fig. 3) is 
caused by an ordered third-layer C configuration; this 
is best seen in Fig. 3(c). This idea was implicitly sug- 

10 1 . From our 



gested by previous experimental work 
first-principles calculations, we found that the geometri- 
cal features of the Si-Si dimers in all 12 cells are essen- 
tially the same. Surface C is invisible, or shows up as a 
missing Si atom, in STM experiments |ll|. This implies 
that one cannot distinguish, using standard microscopy, 
between pure Si(100) or Si(100) with vacancies, and a 
configuration with C in the third layer; the only criterion 
is the long-range order and the change of the reconstruc- 
tion. 
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FIG. 3: Density plot of the C distribution on the surface at 
T=850K; (a) Mc = -6.1, (b) fi c = -6.4eV, (c) ^ c = -6.8 
eV for (c). White corresponds to 4 C atoms, black to C 
atoms. Each panel represents an area of 38x38 nm 2 and is \ 
of our simulation cell. 
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